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I. ABSTRACT 

In this report, Mathematical model for generalized nonlinear three dimensional wave breaking equations was de¬ 
veloped analytically using fully nonlinear extended Boussinesq equations to encompass rotational dynamics in wave 
breaking zone. The three dimensional equations for vorticity distributions are developed from Reynold based stress 
equations. Vorticity transport equations are also developed for wave breaking zone. This equations are basic model 
tools for numerical simulation of surf zone to explain wave breaking phenomena. The model reproduces most of the 
dynamics in the surf zone. Non linearity for wave height predictions is also shown close to the breaking both in 
shoaling as well as surf zone. 

Keyword Wave breaking, Boussinesq equation, shallow water, surf zone. PACS : 47.32-y. 


II. INTRODUCTION 

Wave breaking is one of the most complex phenomena that occurs in the near shore region. During propagation of 
wave from deep to shallow water, the wave field is transformed due to shoaling. Close to the shoreline, they become 
unstable and break. In the process of breaking, energy is redistributed from fairly organized wave motion to small 
scale turbulence, large scale currents and waves. 

Classical Boussinesq theory provides a set of evolution equations for sueface water waves in the combined limit of 

weak nonlinearity (characterized by (5 < I) and weak dispersion (/i < 1) with the raio — = 0(1). the parameters 

T 

represent a wave height to water deth ratio and a water depth to wavelength ratio, respectively. 

It has been shown by numerous researchers that Boussinesq-type equations for varying water depth can describe 
nonlinear transformation in the shoaling region quite well. In the last couple of decades, a lot of research effort has 
gone into improving the predictive c^ability of these equations in the intermediate water-depth and close to the surf 
zone (see e.g. Nwogu @ , Madsen Q, Wei @). It was established that to extend the validity of these equations to 
the deep water, higher order dispersive terms will have to be retained. To improve the predictive capability close to 
wave breaking, the higher order nonlinear terms are very important to include in the equation. However, all these 
models use additional terms that is artificially added to the momentum equation, which would then reproduce the 
main characteristic of a breaking wave, i.e. the reduction in wave height to describe wave breaking phenomena. Wave 
breaking in the software FUNWAVE (FUNWAVE is based on the model described by Nwogu 0) is modeled by 
introducing momentum mixing term developed by Kennedy et al Q. 

Starting with the works of Nwogu @ and Madsen Q, most progress have been done to explain wave breaking 
phenomena. Shen addresses the problem with partially rotational flow. Vorticity dynamics and formation into 
the fluid is very important in the wave breaking as well surf zone. To address this problem, Veeramony & Svendsen 
m derived breaking terms in Boussinesq equation assuming flow as a two-dimensional rotational flow. Here, the 
breaking process is modeled by assuming that vorticity is generated in the roller region of the breaking wave and 
solving vorticity transport equation to determine the distribution of the vorticity. This naturally introduces additional 
terms in the momentum equation which causes wave height reduction as well as changes in the velocity field. However, 
since this model is based on stream function formulation, it cannot be trivially extended to three-dimensional flow. 
The phenomena of wave breaking in Boussinesq equations are being modeled usiM quite few techniques which can 
preserve the wave shape as well as include energy dissipation mechanism. Shen |6| developed a generalized form of 
Bousinesq equation in three dimensional by introducing vertical flow field with arbitrary vorticity distribution up to 
0{pi?). But he did not describe momentum transport equation with full description of rotational flow. Recently, Zou et 
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al @ addressed the problem by including the higher order terms in Boussinesq equation for two dimensinal flow. This 
model solves vorticity distribution based on the parametric form taken form surface roller experimental data. In this 
paper, we developed a general analytic form for breaking term for fully nonlinear set of Boussinesq equations for three 
dimensional vertical flow field near surf zone region. Derivation of breaking term from Reynold stress based vorticity 
transport equation was also developed to describe rotational field as a complete model for thorough understanding of 
wave breaking phenomena. 

The term nonlinear indicates that no truncation based on powers of 6 is employed. The paper is organized as follows: 
Section 2 discusses the basic governing equations for continuity and momentum with boundary conditions. Section 
3 describes the equation for horizontal and vertical velocity distribution for potential and rotational components. In 
section 4, the breaking term is derived for velocity transport equation for fully nonlinear case and solved analytic form 
vorticity transport equation from Fourier series expansion. Reynold stresses are analyzed in the breaking region. In 
last section, results were discussed with conclusion. A parametric analysis of the role eddy viscosity profile is being 
discussed in detail. 


III. BASIC EQUATIONS 


We consider a three-dimensional wave field with free surface rj(x, y, t) propagating over a variable water depth 
h(x,y). As we are primarily concerned with wave breaking, we only consider here wave propagation in shallow 
water. Wave in this region can be characterized by two non-dimensional parameters 6 = a/h and y = h/l where a is 
the characteristic wave amplitude and I the characteristic wave length. The parameter /r is a measure of frequency 
dispersion and 6 that of the non linearity of the wave. In this study, since we are only considering shallow water 
waves, we only have to consider weakly dispersive waves (up to but have to retain all nonlinear terms. In this 

paper, the variables are non-dimensionalized using following scaling: 


u = 


X = x/l,y = y/l^z = z/h^ t = 

u,v = V^W = w 


( 1 ) 


where the ~ represents the dimensional variables, g is the acceleration due to gravity, u and v are the horizontal 
components of the velocity in the x and y directions respectively, w is the vertical velocity. We start with the Eulerian 
equations of continuity and momentum in nondimensionalized form for velocity field u = ( m , v, w) as: 


du du du du dp 
— + Su— + 6v— -I- Sw— + — 
ot ox ay oz ox 


dv 

m' 


Sv 


dv 

dx 


5v 


dv 

dy 


Sw 


dv 

dz 


dp 

dy 


5p^—+ S^p^u— + S^p^v— + S^p^w— -b (5^ -f 1 
ot ox oy oz oz 


Since the fluid flow is rotational, we also have three dimensional vorticity field s 


as 


0 

0 ( 2 ) 

0 


{sx,Sy,Sz) in the fluid defined 


V X u = s (3) 

where V = {d/dx,d/dy,d/dz). The continuity equation then becomes, 

dw , , 

y-u+— = Q (4) 

Here - u = {du/dx, dv/dy). The above equations satisfy two boundary conditions for velocity at bottom and at free 
surface. At the free surface z = r](x,y,t), since particles are free to move with fluid velocity, the kinematic boundary 
condition is 


dp 

Wr,=Uy\7p+ — 


( 5 ) 


and at bottom z = —h{x, y) 


Wb = -Ub ■ \jh 


( 6 ) 
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where = {uri,Vrj) is two component horizontal surface velocity. = {h^,hy) refer to horizontal 

derivative with respect to x and y in all subsequent calculations. The horizontal component for vorticity field s = 
(sy, —Sx) can be described as, 

du „ 

— -fi\w = s (7) 

with u = (u, v) as two component horizontal field whereas vertical component of vorticity expressed as 


du dv 
dy dx 


( 8 ) 


This is straightforward calculation from equation (6) and (8) which is the beginning equation in three dimensional 
vorticity field formulation. 


w + 


d\ 


dz"^ 


= -V ■ s = Su 


(9) 


w represents the vertical velocity of the flow. In the above equation, once w solved, horizontal component of velocity 
u,v can be solved from vorticity relation. In weakly hydro static case ( 0 < ^ 1 ), solution is typically obtained 

from iterative perturbation procedure with successive correction term up to 

In case of breaking waves where vorticity is very strong, so {du/dz ^ 0(1)). We assume solution as, u = Uo + 
and w = wq + for horizontal and vertical velocity component. Under this assumption, 

Poisson equation becomes 


5^2 


( 10 ) 


d'^wi 


d^wo d^wo 
dx'^ dy"^ 


( 11 ) 


wq, w\ can be calculated from bottom boundary conditions using equation (7) separately where the boundary condi¬ 
tions are, 


Who = -Uho ■ \/h 
Wbi -b Ubi ■ Vl*- = 0 


( 12 ) 


at bottom boundary z = -h . 

Since at any other depth z = Zp, w is constrained by continuity equation only, so the equation follows 


dw du 

■^]zr = “V ■ Um]zr + ■ V^r]zj. (13) 

where Um is velocity at any arbitrary depth Zr ■ In Boussinesq equation, one may take depth average or any inter¬ 
mediate velocity for horizontal velocity between bottom and free surface as reference velocity. In the wave breaking 
zone where the vorticity is developed non uniformly, the equations become simpler with the choice of depth average 
velocity which includes contribution from surface vorticity gradient. We assume solution for velocity comes also from 
rotational contribution due to vorticity at the wave surface. So the velocity has both potential as well as rotational 
component, u = Up + Ur , w = Wp + 

We solve wo,wi and uq, ui at any depth z^ 


dwQ 

dz 


-V ■ {Um - ZrS]zr) + ' « 


(14) 


= [VW^o]z^ ■ VZr 


(15) 


duo dui 


Vwo 


with boundary condition [uqJzi. = u^ and [uijzr = 0 

Equations (4) - (16) form basic shallow water Boussinesq equations. 


( 16 ) 
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IV. EQUATION FOR HORIZONTAL VELOCITY 


In the surf zone, vorticity grows very strongly as a non uniform function over depth. Following Shen [2000], we 
define reference velocity as ft = u + Art — r]Sri in terms of depth average velocity u and magnitude of vorticity at free 
surface Sr/ with the assumption of V ' ^ 0. we set here = rj ss linear calibration for Zj. = rij] h) — h does not 

hold here in presence of nonuniform velocity as wave dispersion properties change both spatially and temporally with 
vorticity. And boundary condition can be set as 


dw 

dz 


— Ir, = V M + ?7(V ■ 


Integrating equation (9) from bottom to surface and applying boundary condition to (16) we get wq as , 

Wo = Wbo - (-V -u + VV ■ - S^o 


(17) 

(18) 


where 


SwO 



V ■s)dzdz 


is the vertical velocity distribution generated by horizontal divergence of vorticity added to surface velocity. Now, 
once Wq is calculated, Ui can be calculated from eqn (15) with surface boundary condition {mq},, = and {ui},; = 0 
Finally, we calculate horizontal velocity as 


r 

u{z) = Ur) — sdz + /r^(S'wi — •S'wi) 

J Z 

+ - ^Z^)V(V -U-VV-Sr,) 

- i7z) [v((m + VSr,) ■ S/h) + (v U - V 'Sr,) V ^] + O(m^) 
which on averaging over depth yields, 

- 2 

U = Ur^ - Au + V [v ■ W - ’IV • St?]-[V^ + VSr,] ' 

-(v -u-rjx/ ■ Sr,)\/h + 0(/) 


(19) 


( 20 ) 


Au = A u{z)dz is the average surface velocity contribution due to vorticity and it is significant for suspended 

Iiri 

sediment particles in the flow. The term Au(z) = sdz is the change due to depth variation of vorticity s. The 
total water depth and surface elevation are taken as = z + h and = r] + h. The contribution for velocity 
has and rotational component apart from potential due to vorticity generation. 

After we redefine Hr, = d and z = H^/Hr,, we express potential and rotational component up to order as, 


^2 2 

Up{z) =Up + y(- - z^)d‘^ V (v ■ Up) 
- z)d [v(Mp • Vh) + (v ■ Up) V h] 

Uriz) = Ur - Au{z) + rjSp + - Srul) 

z^)d^VV iS7 ■ Sr,) 
~ [(V ■Sr,)s/h- \j{Sr, ■ VM] 


Similar expressions for vertical velocity are 


Wp{z) = -{h + z) V •'Wp(z) 

lA 1 

= - V • [(^ + z)up] - y((- - z^) V + z)(v(v • Mp)] 
-d^{\ -z)[y-{h + z){s/{up • syh) + (v ■ Up) v h)] 


( 21 ) 


(22) 
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Wriz) = - \/ -{h + z)Ur — V ' [(^ + ^)VS7j] 

r 1 

-y V-ih + z){--z^)d^\7ivV-s-n) 

+Ai^ V -[(^ - z)d [v(??Sr,. V ^) - (V • ’nsv) V h] 


(23) 


V. BREAKING MODEL [OUR CASE: FULLY NONLINEAR] 


The most obvious approach to solve Boussinesq equation is to drop notion of pursuing an expansion in powers of 5 
and instead use weakly dispersive expression for cj) or horizontal velocity in the form of power series in fj? to evaluate 
complete surface boundary condition. We refer to this procedure as fullynonlinear in the sense that all of the availble 
information on velocities is used to evaluate boundary conditions. Conventional time dependent Boussinesq equations 
for surface wave height and consequent breaking term calculation are very straight forward and have benn calculated 
in case of irrotational waves. Here we take up fully nonlinear calculation as vorticity becomes a large fraction of water 
depth in the surf zone or shoaling waves. While developing Boussinesq equations for horizontal momentum, we retain 
up to order 0(:5^) and 0{S^^) in our fully nonlinear calculation. Fully nonlinear Boussinesq equations for long wave 
have been derived by Mei [1983] for flat bottom and by Wei et al Q for variable bottom surface in case of irrotational 
wave. Shen Q addressed problems in developing generalized three dimensional irrotational propagating wave field to 
include rotational motion in general but did not describe the vorticity breaking terms. For horizontal propagation of 
waves, the three dimensional problem can be reduced in terms of two horizontal velocity by integrating over depth 
and retaining up to order 0(6'^) and 0{Sfj?) As horizontal velocity is governed by momentum equation at the surface 
V by, 


Du^ ,Dw, 

-k = (ywI’J + 


Dt' 


Dt 


(24) 


In the surf zone region of sloping beach, waves break due to high vorticity and the breaking of wave later being 
converted to turbulence. So horizontal variation of water depth h{x, y) must be considered in this case. We express 
surface propagation equation in terms of average velocity description and total time derivative of horizontal momentum 
can be written as, 


Du 




du 


+ Ur^ ■ (VW)!,, 


where surface velocity is given by. 


2 2 
_ 2 ^ 

Ur,=u + riSn- —d SJ {y ■ u - 'qsj ■ Sr,) + y^^iV (« - Au|_h + rysr?) ' V^} 


We consider \/Hr^ = yq + yh for wavy bottom 


(25) 


(26) 


Du du dsrj 


(27) 


Du du dSri 2 r 



S + i ■ V} + O)/) 

This long wave momentum equation upon simplification over flat bottom case can be compared to the one derived 
by Shen [2000] The vertical velocity can be obtained similarly. 


Dw 


Dt'^ 


dw dw 

-^\v + Ur, ■ ywn + W—lr, 


(29) 
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So, we can write the horizontal momentum equation as, 


dU M 9 r / N 

— + ti • VM + V»7 = yc? {V(V • m) • Vu + v(V ■ y) 

du du 

+(n • v(V • m)} - ^ V (V • y) + (V • 

-u-V(V- m)}V?? 


(30) 


u is defined in previous section. In contrast to the result by Shen additional contribution factor here arises from 
vorticity variation which is significant for surf zone wave. Wei et al Q also calculated breaking term for irrotational 
long wave momentum equation over a variable bottom wave. The intermediate depth velocity Za is being used there 
proportional to h instead of depth average velocity used here which may not be valid inside the fluid. The use of in 
our approach avoids this difficulty. Finally we try to generalize equation by solving vorticity from vorticity transport 
equation in next section. 


VI. VORTICITY TRANSPORT EQUATION IN BREAKING ZONE 

Madsen and Svendsen Q used a cubic vertical distribution of rotational velocity based on roller jump data which 
can not considered in three dimension case as it is not guaranteed to bring accuracy in the simulation. So we try to 
solve vorticity function from Reynold stress based equation. 


du 1 

a + (“^v)« = --w 


(31) 


Taking the curl on both sides and use vorticity function s = x u we get. 


ds 

y - (s • V)s + (w • v)s = i^V^s 


(32) 


(s • v)® is called ’’vorticity stretching” factor and is the gradient in vorticity value. This term leads to change of 
rotation of material particles present in the flow to the beach. Contribution of this term can not be incorporated from 
two dimension roller jump data. 

We generalize the equation in three dimension as 


ds 


ds 

dx 


ds 

dy 


ds 
’ dz 


du 

dv 

^ dw 

Ss- - Ss — 

— Ss-— 

dx 

dy 

dz 

di^s 

^d^s 

d'^s' 


dx^ 


(33) 


' dy^ dz^_ 

After changing the variable from (x,y,z,t) to wave following coordinates {x,y,a,t)^ we write the vorticity equation 


as 


ds 


dt {h + Sr]) dt da 


Sa dr] ds x . / x ^ 

+ ou(v ■ s) - os(v ■ u) - 


dw ds 
s— - w 


-S^ 


[h + S-q) [ da da 

au , ds „ as ^ ^ du 

(V ■»?) — + 77—TZt(V ■ V)-^ 


s + 


{h + 5q) da (h + 5q) 
1 d^s 


(34) 


{h + Sq) da^ 

The boundary conditions in new coordinate system are, 


+ 0(/i ) + 0{h^) + 0{hy) 


s{a = l,t) = s{x,y,t); 
s{a = 0,t) = 0; 
s{a, t = 0) = 0 


(35) 
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After we redefine s = O + auJs-, which transforms the equation, 


with new boundary, 


dQ dujs ^ O' dr] drjdfl 

dt ^ dt {h + Si]) ^ dt dt do 
So diSt So^ 

Su d^l 6w 
(h + Srj) do~^ {h + Srj) ® 


Sn dw SoLUg dw 
{h + Srj) do {h + Srj) do 


—(50(v ■ u) — SooJs(\/ ■ u) + 


s^on 

{h + Srj) 


iv-v) 


du 


S'^ofl du 

{h + St])^"^ ' 


S^oojg du 

^ {h + Sr])^^ 


(36) 


V,{o = 1, t) = 0 

n{o = o,t) = 0 


(37) 


with initial condition VL{o,t = 0) = 0. This additional equation can be solved numerically as done by Briganti et 
al. [l| for the two dimensional case or an analytical solution can be formulated as shown by Veeramony & Svendsen 
m- The analytical solution can be calculated by assuming n = + Su^'^'^ which gives first and second solution as 

0(1) Problem 


duj^^'i ^ _dujg V 92a;(i) 
dt “ h2 Qjp. 


where the solution is 


(38) 


p(i) 

n 


(- 1 )" 


2 5ws 
UTT dt 


(39) 


duj 

assuming — X)n=i°°'^n sin(n7r(T) And to solve assume = YiriG^\in{m:o) which gives zeroth order 

solution as 


= (- 1 )"— f 

UTT Jo dr 


To consider O(^) Problem 


do h do"^ 


— ^(2) 


where 


odrjduj^^^ o'^ dr] o 

-- w V • : 




' U KJ'JJ U /-i \ 

hat aa ft a-> 

udoj^^'^ 


h do 


(40) 


(41) 


(42) 


To solve above equation, assume sin{mTo) where solution becomes 





(43) 























The solution for vorticity s becomes, 


s = aoJs + sin{mTa) + sin{mTa) (45) 

To solve breaking term, we need value of Ws for boundary and eddy viscosity value as input data. The parametric 
form of vertical profile of eddy viscosity is being applied to calculate vorticity function using various experimental 
data. The eddy viscosity distribution over the water column N(z) is assumed such that its maximum value is located 
at the water surface except at the roller where the maximum is located at the lower edge of the roller. The value is 
estimated by a mixing length hypothesis and r't is given by, 

vt{z,x) = iyoth{x)y/gh{x)N{z) (46) 

The best fit for Wg based on experimental data for hydraulic jump Q, 

X - Xt 

, X — Xt^ , “^0 i , , , 

u;g = 15.75(1- ^)(l-e Ir ) (47) 

tr 

where Xt is the position of the toe of the roller and 1^ is the roller length. Figure 1. shows our analytical calculation 
time evolution of vorticity profile for = 0.03. 


Vorticity Time Profiie 



FIG. 1: Figure 1. shows analytical calculation for time evolution of vorticity profile 


VII. CONCLUSION 

Finally we conclude here by developing a most generalized form fully nonlinear Boussinesq equations for wave 
propagation in surf zone region with variable bathymetry. The vorticity distribution was calculated using Vorticity 
Transport Equation(VTE). In the wave breaking zone, vorticity generated by the shear stress of current is very strong, 
we showed in our calculation the contribution to the surface velocity due to vorticity variation which has significant 
contribution in fluid flow. These extra terms in generalized equation complicate the numerical technique as these 
terms are present in the equation in multiple form of equations for vorticity components which has to be solved 
in coupled solution technique. Veeramony [l0| used simplified formulation by taking constant eddy viscosity value 
but this oversimplihed case may bring inaccuracy in calculation. Briganti et.al [l| formulated a numerical technique 
scheme to solve VTE using generalized depth variable eddy viscosity u = v{x, y) in two dimension case. In three 
dimensional formulation, the nonlinear terms in the vorticity transport equation(VTE) will complicate the calculation 
and so proper numerical technique have to be developed. Turbulence develops because of the instability of vortical 
flow, turbulence injection can be explained from vortical function development. 


VIII. LIST OF SYMBOLS 

a - Wave Amplitude; I - Characteristic wave length; s - Vorticity function; g - gravitational Constant; / - bottom 
friction constant 
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